% Brief Introduction:
% This function converts the B-spline coefficient vector to
% PP-representation coefficient matrix by solving the equation of the
% 0,1,..,k order derivatives at first (#breaks-1) break points.
% Usage:
%   [PCoef,BasisMat]    =   BCoef2PCoef(breaks,k,BCoef)
% Inputs:
%   breaks      :   Vector of break points, accendingly sorted
%   k           :   Polynomial order of the spline
%   BCoef       :   (Optional) The B-spline coefficient vector, 
%                   its dimension is (#breaks+k-1) x 1
% OUTPUTS
%   PCoef       :   Polynomial coefficient matrix. 
%                   The dimension is (#breaks-1)x(k+1) and the spline
%                   function in interval i can be represented as:
%                   \sum_{j=1}^{k+1}(x-breaks[i])^{k+1-j}
%   CoefTrMat   :   The transformation matrix from B-spline coefficient 
%                   vector to the stacked Polynomial coefficient vector.
% Notes:            When BCoef is provided in the function input, the
%                   output will be [PCoef,CoefTrMat], otherwise, the output
%                   will be [CoefTrMat].
% Author:       Xing Guo, University of Michigan, xingguo@umich.edu
% Version:      Nov 28, 2017               
% USES: Spline1dBasis

function [varargout]=BCoef2PCoef(breaks,k,BCoef)

% Compute the Basis matrix at the first (#breaks-1) points at order 0~k
UniBreaks       =   unique(breaks);
BreakBasis      =   Spline1dBasis(breaks,k,UniBreaks(1:end-1),(k:-1:0)');

% Construct the coefficient transformation matrix
for i=1:k+1
    BreakBasis{i}   =   BreakBasis{i}/factorial(k+1-i);
end
CoefTrMat       =   cat(1,BreakBasis{:});
if nargin<=2
    varargout{1}    =   CoefTrMat;
    varargout{2}    =   UniBreaks;
else
    % Reshape the stacked polynomial coefficients
    PCoef           =   reshape(CoefTrMat*BCoef,[length(UniBreaks)-1,k+1]);
    varargout{1}    =   PCoef;
    varargout{2}    =   CoefTrMat;
    varargout{3}    =   UniBreaks;
end